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Abstract: Orbit determination (OD) analysis results are presented for the Lunar Reconnaissance 
Orbiter (LRO) using a commercially available Extended Kalman Filter, Analytical Graphics ’ 
Orbit Determination Tool Kit (ODTK). Process noise models for lunar gravity and solar radiation 
pressure ( SRP ) are described and OD results employing the models are presented. Definitive 
accuracy using ODTK meets mission requirements and is better than that achieved using the 
operational LRO OD tool, the Goddard Trajectory Determination System (GTDS). Results 
demonstrate that a Vasicek stochastic model produces better estimates of the coefficient of solar 
radiation pressure than a Gauss-Markov model, and prediction accuracy using a Vasicek model 
meets mission requirements over the analysis span. Modeling the effect of antenna motion on 
range-rate tracking considerably improves residuals and filter-smoother consistency. Inclusion of 
off-axis SRP process noise and generalized process noise improves filter performance for both 
definitive and predicted accuracy. Definitive accuracy from the smoother is better than achieved 
using GTDS and is close to that achieved by precision OD methods used to generate definitive 
science orbits. Use of a multi-plate dynamic spacecraft area model with ODTK’s force model 
plugin capability provides additional improvements in predicted accuracy. 
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1. Introduction 

The Goddard Space Flight Center (GSFC) Flight Dynamics Facility (FDF) has performed orbit 
detennination (OD) and related mission support for the Lunar Reconnaissance Orbiter (LRO) 
since LRO launched in June 2009. The LRO mission achieved an initial elliptical commissioning 
orbit after lunar orbit insertion. A series of descent maneuvers were perfonned to navigate LRO 
into a 50-km circular orbit around the Moon, where it remained from September 2009 until 
December 2011. In December 2011, LRO was returned to an elliptical (40 km x 180 km) orbit, 
similar to its early mission commissioning orbit. 

Since LRO launch, the FDF has employed the Goddard Trajectory Detennination System (GTDS) 
for operational LRO OD. GTDS is a batch least-squares (BLS) estimator, employing high-order 
gravity modeling and a spherical spacecraft model for solar radiation pressure (SRP) as the 
primary forces for routine OD. In support of mission science objectives, the requirements 
established for LRO OD were: 
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• Definitive ephemeris accuracy of 500 meters total position root-mean-squared (RMS) and 
18 meters radial RMS, 

• Predicted orbit accuracy less than 800 meters root-sum-squared (RSS) over an 84-hour 
prediction span. 

Definitive and predictive accuracy requirements for LRO were easily met using GTDS during the 
50-km circular orbit mission phase. However, predictive accuracy has been poorer in the elliptical 
commissioning orbit, both following initial lunar orbit insertion and currently, particularly during 
high beta-angle periods. Previous work [1] has shown that the use of a more detailed spacecraft 
area model, with definitive attitude modeling for the spacecraft and solar array, can greatly 
improve prediction accuracy in the elliptical orbit. 

This paper examines the use of a commercially-available Extended Kalman Filter (EKF), the Orbit 
Detennination Tool Kit (OTDK), for LRO OD. ODTK implements an EKF and fixed interval 
smoother with high-fidelity force modeling for OD around the Earth or other central bodies, and 
supports over 150 different observation types, including the Universal Tracking Data Format 
(UTDF) range and range-rate observations used in this study. Sequential range and total count 
phase measurements from the NASA Deep Space Network (DSN) are available in small quantities 
for LRO and supported by ODTK but were not used in this study. 

An EKF has a few natural advantages over a BLS estimator for LRO OD. For example, an EKF 
delivers faster OD processing, because an EKF only processes the measurements once rather than 
multiple times. This can save a significant amount of processing time when a high-order gravity 
model is used. In addition, the filter can provide a covariance estimate that is timedependent and 
potentially more realistic than is achievable by GTDS because of the inclusion of dynamical 
process noise. Furthermore, ODTK implements, via a plugin programming interface, a more 
capable box and wing spacecraft area model than currently available in GTDS. This offers an 
opportunity for improving the predicted accuracy for daily operations. 

2. Dynamical Process Noise Model 

When used in an optimal manner for operational OD, the Extended Kalman Filter can be described 
as a continuously running recursive machine. Although operational workflows typically involve 
running the OD process at one or more discrete times each day, each run of the EKF starts with 
final state estimate and state estimate error covariance from the prior run. The resulting filter 
estimates are, therefore, equivalent to the case where the filter is run continuously. Starting each 
OD process with the final state and error covariance from the prior run is desirable since it is 
operationally straightforward and allows for current time estimates to reflect the maximum 
possible amount of measurement information. 

The recursion inside the EKF consists of an alternating series of time updates and measurement 
updates. Time updates are perfonned to advance the state and state error covariance estimates to 
the next measurement time. In OD, the time update includes the nonlinear propagation of the orbit 
and stochastic parameters and the linear propagation of the state error covariance. Dynamical 
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process noise is incorporated during the time update and can be described as additional uncertainty 
injected into the state error covariance to account for the inaccuracies in the model of the dynamics 
of the system. The inclusion of additive dynamical process noise in the covariance portion of the 
time update takes the form 


Pk+ 1 \k = tyk,tk+l)Pk\kQr(tk,tk+l) + Q (tk,tk+l) , (1) 

where Pk+ i \k is the state error covariance at time 4+1 including measurement information through 
time tk, Q{tk,tk+\ )is the linear state error transition function and Q ( 4,4+ 1 ) is the additive process 

noise matrix, which accounts for dynamical uncertainty over the interval from time 4 to time 4 + 1 . 
As seen in Eq. 1, the dynamical process noise matrix is added to the propagated state error 
covariance at time 4 + 1 . 

Measurement updates are performed to fold in measurement information, producing an improved 
state estimate at the time of the measurement and a corresponding reduction in the uncertainty of 
the estimate as expressed in tenns of the state error covariance. It is the role of the dynamical 
process noise to ensure that the covariance is realistic at the time of measurement processing. If 
the dynamical process noise is too small or nonexistent, the filter can become smug. A smug filter 
estimate is characterized as having an error covariance that is significantly too small, which leads 
to the excessive rejection of measurement data and eventually to filter divergence. If the dynamical 
process noise is too large, the automatic data editing of the filter may not be able to effectively 
discern good data from bad data, which will result in degraded state estimates and may also lead 
to filter divergence. For the filter to perform optimally, the dynamical process noise must be 
specified to appropriately account for deficiencies in the model of the dynamics. 

While dynamical process noise can be introduced into the time update in many ways, the preferred 
methodology in ODTK is to associate dynamical uncertainty directly with the physical models 
used in the propagation of the spacecraft trajectory. Wright [2] lists this relationship in his 
requirements for optimal OD. When such “physically connected” process noise is used, the filter 
becomes more robust to variations in orbital conditions, provides a more realistic measure of orbit 
uncertainty and requires less operator intervention. Notable sources of dynamical model 
uncertainty for most space missions include the modeling of central body gravity, atmospheric 
drag, and SRP. In the case of LRO, atmospheric drag is not significant leaving the lunar gravity 
model and SRP as the main sources of dynamical model uncertainty. 

The uncertainty in historical lunar gravity field solutions has been a limiting factor in achievable 
orbit accuracy [3]. Large uncertainties in lunar gravity field solutions were reflective of the lack 
of tracking data for satellites on the far side of the moon for use in gravity model development. 
While the uncertainty in lunar gravity has been greatly reduced with the introduction of the Gravity 
Recovery and Interior Laboratory (GRAIL) mission derived gravity solutions [11, 15], it is still 
important to model the remaining uncertainty in the gravitational acceleration to achieve optimal 
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estimation perfonnance. Accordingly, a process noise model based on the GRAIL gravity solution 
has been developed. The gravity process noise model is based on the integration of acceleration 
errors averaged over a sphere where the radius of the sphere is derived from the location of the 
spacecraft [4, 5]. The output of the gravity process noise model is a 6x6 symmetric matrix, which 

is inserted into the appropriate location in Q ( 4,4+ 1 ) , see Eq. 1. 

Q {tk,tk+ 1 ) incorporates the increase in the orbit state uncertainty from the integration of the 
averaged acceleration errors over the time update interval. 

The lunar gravity process noise model is comprised of two main components: modeling of errors 
of commission and modeling of errors of omission. Errors of commission are errors in the values 
of the gravitational coefficients used in the computation of the gravitational acceleration. Errors 
of omission are errors experienced in the computation of the gravitational acceleration from the 
use of a truncated model of the gravity field. A spacecraft in orbit about the Moon can travel at 
very low altitude allowing its motion to be influenced by the structure of the mass distribution in 
the upper layers of the lunar mantle. Accurate modeling of the acceleration from lunar gravitational 
forces therefore requires the evaluation of the spherical harmonic potential function to a high 
degree and order. The intense nature of such computations motivates the truncation of the available 
fields to a lower degree and order as a matter of practicality. The increase in error in the dynamical 
model as a result of truncation is accounted for in the gravity process noise model as errors of 
omission [4]. 

The gravity process noise model implemented in ODTK is driven by a set of precomputed inputs 
generated from the formal error covariance associated with the gravity field solution [6]. The 
generation of the gravity process noise model inputs for the GRAIL derived gravity fields followed 
the same general procedure that has been validated during the generation of multiple gravity 
process noise models for Earth gravity models [7], but required minor updates in the 
implementation to account for the much higher degree and order of coefficients in the GRAIL 
gravitational field solutions. At the heart of the gravity process noise computation is the evaluation 
of a set of double integrals where the outer integrals are evaluated analytically using an 
approximation enabled by the mean value theorem and the inner integrals are evaluated using pre- 
computed polynomials in the orbit radius, which are specific to the degree of truncation in the 
gravity field [6]. The time required for generation of the inner integral polynomial representations 
grows significantly with increasing degree and order of the gravity field for several reasons. First, 
the integrand requires the computation of a sum of Legendre polynomials to the degree of the 
gravity field. Second, polynomials are generated for each possible user selected degree of 
truncation. Finally, higher degree gravity effects add higher frequency content to the integrand, 
which subsequently requires smaller step sizes in the numerical quadrature to produce an accurate 
result. Fortunately, these burdensome computations are only required in the pre-processing step 
and do not affect the performance of the filter during estimation. 

SRP induced accelerations are analyzed using two different models for comparison purposes and 
to detennine the importance of detailed spacecraft modeling in the prediction of the orbit of LRO. 
The first SRP acceleration model assumes a uniform spherical spacecraft body augmented with a 
stochastic parameter to allow for variations in the magnitude of the acceleration along the Sun 
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line. The stochastic scale factor of acceleration for SRP accommodates the injection of dynamical 
process noise along the Sun line. This “along-axis” process noise is incorporated by setting the 

diagonal element of Q ( 4,4+ 1 ) associated with the SRP scale factor to a value computed based on 
the selected stochastic sequence for the given time update interval. In addition, ODTK provides 
the option to include fractional white noise in directions orthogonal to the Sun line to account for 
the effects of a non-spherical spacecraft body. This “off-axis” SRP process noise is specified by 
the operator as a fraction of the nominal acceleration along the Sun line. These “off-axis” SRP 
acceleration uncertainties are integrated over the time update interval using a localized 
approximation of two body motion to yield a 6x6 symmetric matrix, which is added to the gravity 

process noise in Q ( 4 , 4 + 1 ). We note that the use of white noise in the directions orthogonal to the 
Sun line is appropriate for cases where a changing attitude profile allows for such off-axis 
accelerations to average out, but is less so for cases where off-axis accelerations are significant 
and slow to change in direction. 

The second model for SRP is a more detailed physical model of the spacecraft implemented using 
the ODTK SRP plugin interface. The particular solar pressure plugin used for LRO implements a 
box-wing model of the spacecraft where each surface has an associated area, specular and diffuse 
reflectivities, and orientation in the body frame of the satellite or a frame aligned to the Sun vector. 
This model can directly compute SRP-induced accelerations in directions on and off the Sun line. 
Furthermore, model parameters can be selected for estimation and be associated with stochastic 
descriptions to allow the addition of dynamical process noise related to uncertainties in the 
reflective characteristics of modeled spacecraft components. Dynamical process noise for the box- 
wing model is incorporated by setting the diagonal elements of Q (4,4+1 ) associated with the 
estimated characteristics of the model to values defined by the selected stochastic sequences for 
the given time-update interval. A feature of dynamical process noise that is physically connected 
to the SRP acceleration is that no process noise is added when the satellite is in eclipse. 

3. Analysis Time Span and Methodology 

To facilitate comparisons between GTDS and ODTK performance, the analysis time span chosen, 
March 1 1, 2013 to July 13, 2013, is the same as used in the previous analysis on LRO OD accuracy 
using GTDS [1]. This time period spans low to high beta angle orbit conditions and includes a 
full-Sun orbit period, five momentum unloads, and one station-keeping maneuver. 

Over this analysis time period, LRO tracking consisted primarily of S-band range and range-rate 
measurements from a NASA tracking station at White Sands, New Mexico (WS1S), and from 
Universal Space Network (USN) stations in Perth, Australia; Hawaii, USA; Weilheim, Germany; 
and Kiruna, Sweden. During this period, the NASA Deep Space Network (DSN) provided 
occasional tracking passes during momentum unloads and station keeping maneuvers, but this 
constitutes a small fraction of tracking data. Therefore, DSN data was not included in the analysis. 

A number of relevant LRO tracking data issues have been described previously [1,8] and for 
brevity are summarized here. They are: 
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• An approximate range-rate bias of -1.0 cm/sec on all USN stations, 

• An approximate measurement time-tag bias of +6 milliseconds on WS1S range tracking, 

• An approximate measurement time-tag bias of -2 millisecond on USN range tracking, 

• Various differential range measurement calibration biases among all stations. 

The analysis automation was configured to run the filter as it would be in an operational daily use 
scenario. After initialization, the filter restarts each day in “warm start” mode from a restart record 
at 12:00 UTC and processes 24 hours of tracking data from the prior day at 12:00 UTC to the 
current day at 12:00 UTC. At the conclusion of each daily filter run, the filter generates a 5 day 
prediction using the filter ending state estimate. Following the filter execution, the fixedinterval 
smoother is run backward in time for a total of 96 hours, starting at the filter end time and ending 
on the 4 th prior day at 12:00 UTC. A 4-day definitive ephemeris is generated by the smoother 
process. No prediction is generated from the smoother. Because the smoother start state and filter 
end state are identical, any prediction from the smoother would be identical to that from the filter. 
This processing workflow matches what would be required for operational use of the filter for 
daily LRO operations with OD beginning at around 7 a.m. EST each day. 

Once the filter and smoother processing are complete, a separate utility runs to perform the 
McReynolds’ filter-smoother consistency test [9]. The McReynolds’ filter-smoother consistency 
test computes a unitless metric at each point of overlap between the filter and smoother that is the 
ratio of the difference between the filter and smoother state estimates to the difference in filter and 
smoother formal covariance. For each point at which the metric is less than 3, the filter and 
smoother state estimates are demonstrated as consistent with their associated covariances and the 
test is passed. Formally, 99% of all tested points are expected to pass the test, however, from a 
practical standpoint passing the test for the vast majority of the time based on visual inspection is 
considered sufficient. On the other hand, persistent excesses of the test metric are an indication 
that some observable effect is ignored, mis-modeled, or poorly tuned. The filtersmoother 
consistency test may be applied to any estimated parameters, including biases and the components 
of the SRP model. 

In addition to the filter-smoother consistency test, a suite of other quality assurance (QA) reports 
and graphs were generated after each analysis series run. These included graphs of the filter 
covariance, values of any estimated parameters, and residual ratios or scaled residuals, which are 
the raw observation residuals divided by the root sum of the covariance and measurement noise. 
Residual ratio plots are very useful for ensuring that the filter remains in a converged state, and 
also provide insight into the magnitude of the actual versus modeled measurement noise. If the 
residual ratios are excessively narrowed around zero, this indicates that the either the measurement 
noise or (less often) the process noise is too large. 

The accuracy of the ephemeris, both definitive and predictive, generated by the filter and smoother 
was assessed using the filter and smoother fonnal covariance and by comparison to precision orbit 
ephemeris files generated by the LRO Lunar Orbiter Laser Altimeter (LOLA) team using the Orbit 
Determination and Geodetic Parameter Estimation (GEODYN II) Program [10]. The comparisons 
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between ODTK and GEODYN ephemerides are described in more detail in Section 4.6. Predictive 
or definitive spans crossing maneuvers or momentum unloads were excluded from the statistics. 

For this study, the quality of each tuning scenario was evaluated primarily based upon its 
performance against the mission requirements, especially the predicted accuracy requirement. 
Nearly all tunings achieved definitive accuracy results well below the requirements, with radial 
RMS definitive accuracy (as measured by comparison to the GEODYN solutions) typically close 
to 1 meter, and RMS total definitive accuracy, similarly measured, between 7 and 20 meters. As 
one would expect a filter with dense tracking to easily perform well for definitive accuracy, 
predicted accuracy - and in particular the number of failures of the 800-meter 3-day predicted 
accuracy requirement - was a critical metric for evaluating the quality of each tuning scenario. 
Nevertheless, in addition to the definitive and predicted accuracy, the value and behavior of 
estimated parameters, in particular SRP and the filter-smoother consistency of both the orbit and 
estimated parameters, were also used as discriminating criteria in evaluating filter tuning. 

4. Filter Tuning 

The major forces affecting LRO OD are lunar gravity and SRP. Second order forces such as lunar 
albedo, thermal radiation pressure, and Earth oblateness were ignored in this study. The goal of 
filter tuning is to develop a filter configuration that runs autonomously or with minimal 
intervention, achieves the best definitive and predictive orbit accuracy possible, and generates a 
realistic covariance. To achieve this goal, the orbit analyst must study and consider all of the 
relevant forces and sources of error and model them more carefully than is generally required for 
operational BLS OD. The filter tuning process requires the analyst to address both the nominally 
modeled accelerations of these forces and their presumed errors. In addition, the filter analyst must 
establish realistic models of all included orbit measurement types, taking into account their white 
noise and potential measurement and time-tag biases. ODTK offers a variety of stochastic models 
that may be applied to estimation of bias and force model parameters. In most cases operator 
experience or empirical testing must be used to select a preferred model and which parameters 
should be estimated versus being applied. More than 70 runs were perfonned over the entire 
analysis span to study tuning options in a parametric fashion. The following sections describe the 
results from these runs and the motivations for the selected modeling parameters of the major 
forces and measurement types. 

4.1. Gravity Model Selection 

The GL0660B lunar gravity model [11] was chosen for this study. Preliminary ODTK runs were 
conducted to determine an appropriate truncation of the model which would provide a manageable 
balance of OD accuracy and computation time. A 100x100 truncation performed poorly, 
producing 43 failures of the 800-meter predicted accuracy requirement. Higher order truncations 
performed much better, with the 150x150 truncation giving 15 predicted accuracy failures and 
200x200 producing 12 failures. Both the 150x150 and 200x200 truncations had better definitive 
accuracy than the 100x100 case, with the 200x200 truncation performing about 1 meter better in 
mean total position accuracy than the 150x150 case. As a 200x200 truncation is currently used in 
operational GTDS OD for LRO and the processing time of about 9 minutes to filter 24 hours was 
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considered manageable, the 200x200 truncation was selected for the ODTK analysis runs. 
Truncations larger than 200x200 are currently impractical for rapid daily OD using GTDS but may 
possibly be practical for use with a filter due to the filter’s need to only iterate once on the 
measurement data. 

4.2. Solar Radiation Pressure Modeling 

SRP is the largest non-gravitational perturbation affecting the LRO orbit and inadequate modeling 
of SRP is the primary cause of large prediction errors for LRO, particularly during high-beta angle 
periods [1]. ODTK supports estimation of SRP using a random walk, GaussMarkov, or Vasicek 
stochastic model. The Vasicek model is a hybrid model that acts like a superposition of a random 
walk and Gauss-Markov model, facilitating separation of long-term and short-term effects [12]. 
For all models the user may also specify, as a fraction of the nominal acceleration from SRP, 
supplemental off-axis SRP process noise both perpendicular to the Sun line in the ecliptic plane 
and nonnal to the ecliptic plane, to account for SRP force model errors caused by using a coarse 
area model. 

Analysis runs were perfonned using both the Gauss-Markov and Vasicek models with various 
levels of additional SRP process noise. Prediction accuracy using the Gauss-Markov SRP model 
was found to be sensitive to the level of additional off-axis SRP process noise. As the process 
noise was increased, the prediction accuracy improved and the number of predicted accuracy 
requirement failures decreased. The best runs in this series used SRP process noise as high as 
200% of the nominal acceleration. Note that the off-axis model assumes that acceleration errors 
can be modeled as white noise. The large percentage of the nominal SRP acceleration used within 
the model is indicative of a violation of this assumption and that actual off-axis acceleration errors 
are likely strongly correlated over significant time spans. 

The Vasicek model was preferable for SRP modeling. While comparable prediction accuracy was 
achievable for various tunings of the Gauss-Markov and Vasicek models, for many cases 
estimation of the correction to the nominal SRP using a Gauss-Markov model showed large 
variations without any apparent converging behavior. In contrast, estimation of the SRP correction 
using the Vasicek model was less variable and exhibited a converging behavior to a stable long- 
term mean value, with some short-tenn variations still apparent. Even when using the Vasicek 
model, off-axis SRP process noise was still required. A run utilizing the Vasicek model and no 
off-axis SRP process noise diverged very close to the time LRO entered full-Sun exposure near 
the end of May 2013. Estimation of the correction to the nominal SRP coefficient using the 
Vasicek model consistently showed convergence to a correction of about -0.71, so the a priori 
coefficient of solar radiation pressure was adjusted for later runs from 1.67 to 0.96. We were able 
to confirm that this value is consistent with the values of coefficient of SRP obtained by the 
GEODYN precision OD team, so this gives additional confidence to the use of the Vasicek model. 
The following figures show examples illustrating estimation of the correction to the coefficient of 
SRP using a Gauss-Markov (Fig. 1) and Vasicek (Fig. 2) model. 
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Figure 1. SRP correction estimate using a Gauss-Markov model 
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Figure 2. SRP correction estimate using a Vasicek model 


It is interesting to note that the filter enables estimation of the coefficient of SRP. Prior attempts 
to estimate SRP with GTDS were never successful, typically yielding inconsistent or unrealistic 
values. GTDS does not account for process noise, therefore gravity model error easily aliases into 
the estimate of the coefficient of SRP. It is conjectured that the estimation of corrections to SRP 
requires a long time period of observation over which the unique dynamical signature can be 
isolated to provide observability. In such a case, the inclusion of off-axis SRP process noise 
provides additional freedom for the estimator to allocate state corrections into the orbit states and 
avoid aliasing localized errors into the SRP state. From Fig. 2, the apparent time to convergence 
for the estimated correction to the coefficient of SRP is about seven weeks, much longer than the 
GTDS BLS fit span of only 1.5 days. However, it is noted that the filter and smoother estimates 
of SRP state may still not reflect physical SRP force model changes. Over the time span, a number 
of discrete changes to the LRO solar array configuration occur at known epochs, discussed in 
greater detail later. Neither the Gauss-Markov nor Vasicek estimates provide clear indications of 
these changes. 

4.3. Range-Rate Data Modeling 

Initial tuning efforts focused on using only the range-rate tracking data in the solution and excluded 
the range data. Prior LRO OD mission experience using GTDS has shown that rangerate only 
solutions are adequate to meet mission requirements. In addition, although the LRO range-rate 
data does have measurement biases, excluding range data simplifies the tuning process by the 
eliminating uncertainties from transponder delay and station range calibration biases, and avoids 
the issue of the anomalous time-tag biases associated with the LRO range tracking. 

Mission requirements specify that range-rate tracking provided by the USN should have an 
accuracy of 3 mm/sec and that provided by the WS1S should have an accuracy of 1 mm/sec, 
1 sigma. For sake of simplicity, all stations were initially configured with 1 mm/sec of 
measurement white noise. With this setting it was quickly apparent that, at least over the analysis 
span, range-rate tracking from the USN Kiruna stations was noisier than the other stations. As a 
result, the range-rate noise for both Kiruna tracking stations was adjusted to 5 mm/sec. 

As previously mentioned, the USN range -rate tracking is known to have an approximate -1.0 
cm/sec bias for all USN stations. Test cases were run that either applied this bias or estimated it, 
as well as the approximately 0.0 cm/sec White Sands range-rate bias, using a Gauss-Markov 
stochastic model. As the bias was well known from the GTDS solutions, an initial uncertainty of 
1 mm/sec was assumed and test cases were run to vary the half-life of the estimated bias. The 
predicted accuracy improved and the number of prediction accuracy requirement failures 
decreased steadily as the half-life increased. Employing a long half-life with small uncertainty 
approximates applying a bias, so it was decided to apply the nominal mean values of range-rate 
biases and not to estimate them at all. 
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Although nominally orbiting the Moon in a nadir-pointing attitude, LRO performs frequent slews 
(about 2 to 8 per week, but sometimes more) to off-nadir attitudes for science observations or 
instrument calibration. These slews were evident as numerous large range-rate residual excursions 
and were easily observable in the scaled residual plots. Aside from these slews, a close 
examination of the range-rate residuals revealed that most passes had a lower-order parabolic or 
sinusoidal signature (Fig. 3). This signature is induced by the motion of the LRO antenna relative 
to the line of sight during the pass. Most LRO tracking is performed using the high-gain antenna 
(FIG A) which is mounted atop a mast extending approximately 2.7 meters from the spacecraft 
center-of-mass. LRO has two additional omni-directional antennas, mounted on the +X and -X 
spacecraft axes, but these antennas are only typically used during certain momentum unloads or 
spacecraft anomalies. 



SCN.USHS Doppler Meas Residuals SCN.USPS Doppler Meas Residuals 
SCN.WS1 S Doppler Meas Residuals 3-Sigma 


SCN.KU2S Doppler Meas Residuals 
-3-Sigma 


A 


Figure 3. Range-rate scaled residuals without HGA motion modeling 

ODTK provides the capability of modeling range-rate effects from antenna motion relative to the 
spacecraft center of mass, a capability not available in GTDS. Because the FDF receives both 
definitive and predicted spacecraft attitude history files, we were able to execute analysis runs 
using the antenna center-of-mass offset and definitive spacecraft attitude history in the range-rate 
measurement modeling. Modeling the HGA motion using attitude data cleaned up the range-rate 
residuals considerably (Fig. 4), yielding residual distributions closer to white noise. A few residual 
excursions remain, particularly during targeted attitude slews, such as near 17:00 UTC on March 
12 in Fig. 4. Most dramatic, however, was the improvement in filter-smoother consistency (Figs. 
5 and 6). Once HGA motion is accounted for, only a couple brief large excursions still remain. 
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This result illustrates the power of the filter-smoother consistency test in the indication of 
modeling deficiencies. 
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Figure 4. Range-rate scaled residuals including HGA motion modeling 
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Figure 5. Filter-smoother consistency without HGA motion modeling 
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Figure 6. Filter-smoother consistency including HGA motion modeling 


4.4. Range Data Modeling 

Inclusion of range tracking in the filter solutions expands the complexity of the tuning process, as 
it requires the analyst to address the modeling of range measurement biases, LRO’s transponder 
delay, and tropospheric delay corrections to the range observations. These separate effects tend to 
be highly correlated which means that there is no easy choice about whether to apply or estimate 
each. In the case of LRO range tracking, there are also time-tag biases affecting the range data that 
are not present in the range-rate tracking. 

With preprocessing of the tracking data, it is possible to separate the range and range-rate 
measurements, enabling application or estimation of time-tag biases solely to the range 
measurements. When this was done, the best solutions using only range data performed 
comparably to the best range-rate-only and combined range and range-rate cases for definitive 
accuracy, but this range-only case was worse for prediction accuracy, producing two failures of 
the predicted accuracy requirement. The best range-rate only cases had no failures of the predicted 
accuracy requirement. The best runs combining range and range-rate data applied transponder, 
troposphere, and measurement time-tag biases and estimated an individual range bias on each 
tracking station. However, this combined data type run yielded definitive and predictive accuracy 
identical to that obtained using only range-rate data. This lack of improvement makes it generally 
undesirable to take on the added complexities of using the range data in the state estimation. 

4.5. Supplemental Generalized Process Noise 

Detennining an ideal level of process noise is perhaps the most difficult part of the filter tuning 
process for cases where a dominant noise source is not present. ODTK provides physically 
connected process noise to account for errors of commission and omission in the gravity field and 
along-axis and off-axis acceleration errors from SRP as described above. In some cases, however, 
additional process noise is required to account for u nkn own and un-modeled forces or forces for 
which a deterministic model is not known or not included in the dynamical model (such as out 
gassing, thermal re -radiation, and lunar albedo). Quantifying process noise in detail is a daunting 
task. To account for such circumstances, ODTK allows the user to specify supplemental 
generalized white process noise not specifically correlated with any force model. In the case of 
LRO, this additional process noise would primarily account for potential deficits in the presumed 
SRP process noise and gravity model error, as well as soaking up effects of the unmodeled second- 
and higher-order forces. 

It was found that the addition of generalized process noise at the level of 1 mm/sec per orbit, 
equivalent to an acceleration of approximately 1.4xlO" 10 km/sec 2 , on all axes (radial, in-track, and 
cross-track) was beneficial to both predicted accuracy and filter-smoother consistency. For a low 
lunar orbiter, this magnitude of acceleration is on the order of SRP [13]. The inclusion and effect 
of the generalized white process noise is similar to what was described earlier for the offaxis SRP 
process noise where the effect on the state error covariance is much smaller than would be seen 
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by modeling a serially correlated error of a similar magnitude. Further increases in the level of 
generalized process noise qualitatively improved filter-smoother consistency even more, but 
degraded both definitive and predictive accuracy. 

4.6. ODTK Definitive Accuracy 

One natural measure of filter definitive accuracy is the filter fonnal covariance. Fig. 7 shows the 
3-sigma filter covariance for the best-performing tuning scenario. The large, off-scale spikes in 
radial and in-track variance are from modeling of momentum unloads or maneuvers at those times. 
The effects of orbit viewing geometry are apparent in both the in-track and cross-track 
uncertainties, with the in-track uncertainty peaking sharply when the LRO orbit plane is viewed 
face-on. 



Figure 7. Position uncertainty, 3-sigma, for the range-rate-only scenario 

As noted in Section 3, the LRO LOLA science team produces precision orbit ephemerides for 
LRO using the GEODYN program [14]. The latest release of GEODYN-derived ephemerides use 
the GRGM900C gravity model [15] truncated at degree 600. These ephemeris files are described 
as having an orbit accuracy of about 10 meters RMS total position error, with about 1 meter or 
better RMS radial error [16]. The GEODYN GRGM900C solutions were employed as the “truth” 
orbit for comparisons to the ODTK ephemeris over ODTK definitive and predicted spans. 
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Definitive comparisons were perfonned between GEODYN solutions and the ODTK 24-hour 
forward filtered definitive ephemeris, a smoother definitive ephemeris smoothed backwards 48 
hours, and a smoother definitive ephemeris smoothed backwards 96 hours. Comparisons between 
the best-case ODTK runs and GEODYN are shown in Tab. 1. For background on the GTDS 
solutions, please see Ref. [1]. 

Table 1, Definitive Accuracy Results 



Radial 
Mean / RMS 
(meters) 

In-track 
Mean / RMS 
(meters) 

Cross-track 
Mean / RMS 
(meters) 

Total Pos. 
Mean / RMS 
(meters) 

LRO requirement 

None / 18 

None 

None 

None / 500 

GTDS BLS using a spherical 
area model and constrained 
plane 

0/2 

-12/35 

0/29 

30/45 

GTDS BLS using a multiplate 
area model and definitive 
attitude and constrained plane 

1/3 

10/25 

0/16 

21/30 

ODTK filter definitive 

0/10 

4/59 

0/16 

14/62 

ODTK 2-day back-smoothed 
definitive 

0/3 

2/13 

0/6 

7/14 

ODTK 4-day back-smoothed 
definitive 

0/2 

0/10 

0/5 

6/12 


Both smoothed ephemeris files have better accuracy than the filter because the smoother 
dramatically improves the accuracy of the definitive ephemeris during periods immediately after 
momentum unloads and maneuvers, when the filter is still reconverging from coarse delta-V 
modeling. Filter radial and in-track performance is worse than GTDS because of the increased 
error in the filter ephemeris during these reconvergence periods. Steady-state (excluding 
reconverging periods) filter definitive accuracy is better than both GTDS cases, with a mean total 
position error of 9 meters, and RMS total position error of 15 meters. The smoothed ephemeris is 
more accurate than the filter ephemeris, as it removes maneuver modeling errors, but the 
improvement between the 2-day smoothed definitive and 4-day smoothed definitive is small. 

A particular challenge to GTDS OD for LRO is the poor observability of certain components of 
the orbit that occurs periodically as the orbit rotates from a face-on to an edge-on view to the Earth. 
In a face-on view, radial and along-track accuracy are poor. In an edge-on view, crosstrack 
accuracy is poor. With GTDS, cross-track errors can be somewhat ameliorated by application of 
a constrained plane during edge-on viewing geometry [1], However, the filter’s dynamic 
covariance has the advantage of naturally retaining knowledge of the orbit plane estimate from 
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good to poor viewing geometry. RMS definitive cross-track error using GTDS without a plane 
constraint was 60 meters, which was reduced to 29 meters when a constrained plane method was 
applied. In contrast, ODTK’s filter RMS definitive cross-track error was 16 meters, which came 
down to 5 meters with 4-day smoothing. 

5. Filter Performance Using a Multi-Plate Area Model 


As noted in the previous analysis of LRO OD using GTDS [1], a box-wing model performs better 
than a spherical model for LRO prediction, at least in the case of BLS estimation. However, for 
simplicity the initial analysis and tuning effort employed a spherical spacecraft model. 

ODTK supports the ability to model a user-defined box-wing SRP model by the use of a force 
model plugin. ODTK supplies plugin points for arbitrary forces, as well as specific interfaces for 
SRP and drag forces. The SRP plugin point allows a single light reflectance area vector to be 
computed and returned to ODTK in a variety of reference frames, including the spacecraft body 
frame. This area vector is then used to compute the SRP force in the inertial frame by 


F, 




( 2 ) 


where Ois the mean solar flux at 1 astronomical unit (AU), c is the speed of light, R s is the distance 

from the Sun to the spacecraft in AU and is the area vector expressed in the inertial frame. The 
area vector can be thought of as incorporating the effects of area, reflectance and directionality. 

For a simple spherical spacecraft model, ^ would lie along the direction from the spacecraft to the 
Sun and have a magnitude equal to the product of the frontal area and the solar pressure coefficient. 
The plugin also allows the user to add model parameters from the plugin into the state space for 
estimation by the filter and smoother. 

A multi-plate SRP model utilizing this plugin interface has been created by Analytical Graphics, 
Inc., (AGI) and was used for this part of the study. This SRP plugin allows the user to define an 
arbitrary number of flat plates by specification of area (in nr), specular and diffuse reflection 
coefficients, and the unit normal vector components in either the body or Sun frame. The body 
frame is defined by ODTK’s current spacecraft attitude definition, in this case using LRO’s daily 
predicted attitude files supplied by the LRO mission operations center (MOC). The Sun frame is 
defined as the +Z axis aligned to the vector from the spacecraft to the sun and the +Y axis 
constrained to the body +Y axis, with +X completing the orthogonal system. 

For each user-defined plate, the plugin will compute the area vector along with partial derivatives 
of the SRP force with respect to elements of the state space. These calculations are simplified 
greatly by a variety of helper functions in the plugin interface. The area vector for a collection of 
plates representing the spacecraft is simplified from [17] to be 
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( 3 ) 


i = £ [(1 - Pi)s + 2 + pi cos /?,)) n,J 4,c°s /?, 

;=i 

where /V is the number of plates, p, is the fraction of sunlight reflected off the plate, 5, is the 

xv 

fraction of sunlight diffused off of the plate, s is the unit vector from the spacecraft to the Syn, , is 

n - 

the unit vector normal to the plate, A, is the area of the plate, and p, is the angle between , and 5 . 
We note that contributions to the area vector are only included for plates with 0° < (3, < 90° [18]. 
Specular and diffuse reflectivity values for the LRO bus and solar array were available from pre- 
mission measurements. The plugin does not consider self-shadowing of the plates. 

The LRO solar array was used in a number of separate modes during the timeframe of interest. 
When the beta angle, (3/, is low, the solar array tracks the Sun with an offset of 30°. For moderate 
and high beta angles, the solar array is fixed at particular positions. The solar array configurations 
as supplied by the MOC are shown in Tab. 2. 

Table 2. Solar Array Management During the Analysis Span 


Start 

Stop 

Description 

3/11/2013 00:00 

3/12/2013 00:00 

Solar array tracks the Sun, with a 30° 
offset from the Sun vector in the body +Y 
axis 

3/12/2013 00:00 

4/15/2013 13:08 

Solar array tracks the Sun, with a 30° 
offset from the Sun vector in the body -Y 
axis. 

4/15/2013 13:08 

5/21/2013 17:34 

Solar array fixed along the body +Y axis 

5/21/2013 17:34 

7/2/2013 13:42 

Solar array fixed along the body +Y axis, 
with a 30° offset toward the body -X axis 

7/2/2013 13:42 

7/13/2013 00:00 

Solar array fixed along the body +Y axis 


This scheme was confinned in AGI’s Systems Tool Kit (STK) using definitive solar array normal 
vector data provided by the MOC. The definitive data includes numerous brief excursions from 
the nominal behavior listed above that were not modeled. The current plugin assumes a constant 
spacecraft configuration and therefore a separate user input file must be supplied at the time of 
each configuration change. This was done through external ODTK automation scripts. 

The plate model as implemented includes a few known limitations. As previously mentioned, self- 
shadowing (e.g. the solar array shadowing the body) is not modeled in the plugin. However, 
analysis in STK using accurate spacecraft models and STK’s Area Tool indicated there was 
minimal self-shadowing in most cases. Additionally, the HGA boom and dish were not included 
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in the plate model. Similarly, analysis using STK’s Area Tool indicated these errors were typically 
on the order of -1 nr or less. Despite these approximations it was found that use of the multi -plate 
area model plugin eliminated the need for supplemental SRP process noise. 


5.1. ODTK Predictive Accuracy 

Predicted accuracy for the best-case ODTK scenarios, again measured against the GEODYN 
solutions, is shown in Tab. 3. The results reported in Tab. 3 are the mean and standard deviation 
(std. dev.) of the daily maximum 84-hour total position prediction error for each run covering the 
analysis period, along with the number of times the 800-meter prediction accuracy requirement 
was exceeded for each analysis series. 

Table 3. Predicted Accuracy Results 



84-hour Predicted Total 
Position Error 
Mean / Std. Dev. 
(meters) 

LRO requirement 

800 

GTDS using a spherical area model 

325 / 325 8 

and constrained plane 

failures 

GTDS using a multi-plate area 

128/84 

model and definitive SC attitude 

0 failures 

ODTK using a spherical area 

381 /237 0 

model 

failures 

ODTK using a multi-plate dynamic 

243/ 116 0 

area plugin and predicted attitude 

failures 


Both the spherical and multi-plate area plugin models met predicted accuracy requirements during 
the analysis span. Mean predicted accuracy using the plugin dynamic area model was nearly 140 
meters better than that achieved using the spherical area model, but not quite as good as that 
achieved using the GTDS multi-plate area model. One possible reason may be that the GTDS 
predictions used definitive attitude data for the solar array and included HGA modeling, also with 
definitive HGA pointing data. Definitive attitude data is not available in time for daily OD, so the 
results of the best-performing GTDS case are not currently achievable in routine operations. Use 
of predicted attitude data, as in the ODTK case, more accurately represents what can be achieved 
operationally. Omitting the detailed SA behavior and the HGA entirely represent a somewhat less 
complex and easier to configure setup. 
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6. Conclusion 


The results in this paper demonstrate that ODTK’s Extended Kalman Filter can be tuned to yield 
definitive accuracy better than GTDS batch least-squares solutions and predicted accuracy as good 
as or better than the operational GTDS OD. Modeling the effects of HGA antenna motion on 
range-rate measurements in the OD process was vital to achieving good filter-smoother 
consistency when processing range -rate data. As the analysis converged on tunings that gave good 
definitive and predictive accuracy, there was no improvement in the filter-smoother consistency 
until HGA motion was accounted for. This supports the use of the filter- smoother consistency test 
as a sensitive metric for evaluating filter performance. OD definitive accuracy using predicted 
attitude data was equivalent to that using definitive attitude data, which makes HGA antenna 
motion modeling a viable option for operational daily OD. GTDS does not have the capability to 
account for antenna motion in range-rate measurement modeling. Estimation of an SRP 
coefficient, previously not possible with GTDS, was successful using a Vasicek stochastic model 
in ODTK. 

In addition to improved definitive and predictive accuracy over GTDS and the availability of a 
time-dependent covariance, orbit processing using ODTK was faster. In a test of a typical daily 
OD processing scenario, GTDS took 22 minutes to perform state estimation and generation of a 
5 -day predicted ephemeris, requiring 4 iterations on a 1.5 -day measurement arc. ODTK took only 
9 minutes to filter forward 1 day of tracking, generate a 5-day predicted ephemeris, and smooth 
back 4 days. Although GTDS was originally selected as the primary tool for LRO OD based on 
its history as a trusted tool in the FDF and its heritage with the Lunar Prospector mission, ODTK 
is a viable alternative for operational LRO OD and is in a number of ways more advantageous for 
use than GTDS. 
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